PROCEEDINGS 

OF SCIENCE 



o 

(N 

> 
O 



The imbalanced Fermi gas at unitarity 



C/2 

cr 
3 



O. Goulko 

DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 OWA, UK 
E-mail: |p . Goulko@damtp . cam .ac.uk 



M. Wingate 

DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 OWA, UK 
E-mail: \a . Wingate @damtp . cam . ac .~ 



T3 

o 
o 



> 

(N 

m 
o 



Lattice field theory is a useful tool for studying strongly interacting theories in condensed matter 
physics. A prominent example is the unitary Fermi gas: a two-component system of fermions 
interacting with divergent scattering length. With Monte Carlo methods this system can be studied 
from first principles. In the presence of an imbalance (unequal number of particles in the two 
components) a sign problem arises, which makes conventional algorithms inapplicable. We will 
show how to apply reweighting techniques to generalise the recently developed worm algorithm 
to the imbalanced case, and present results for the critical temperature, the energy per particle, 
the chemical potential and the contact density for equal, as well as unequal number of fermions 
in the two spin components. 
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1. Introduction 

With the discovery of the renormalization group came the birth of both lattice gauge theory 
and statistical field theory as we know them today. Nowhere is this kinship more obvious than 
in the study of critical phenomena. For example, the numerical studies of deconfinement in gauge 
theories are identical in most respects to those of demagnetization in spin systems. Now a challenge 
in both fields is to develop more efficient methods for systems with strongly interacting fermions. 
Then we could better understand the phases of gauge theories with matter and of superfluids and 
superconductors. 

Given that the common problem lies in dealing with fermionic degrees of freedom, we are 
motivated to consider the unitary Fermi gas, which is beautiful in its simplicity. This is a system of 
2-component nonrelativistic fermions (typically ultracold 6 Li or 40 K atoms in experiments) inter- 
acting via a short-range potential. The gas is so dilute that only s-wave scattering is relevant, thus it 
is sufficient to treat the interaction as a local 4-fermion coupling. On a spatial lattice, the simplest 
Hamiltonian is that of the attractive Hubbard model. 

The lattice theory describing unitary Fermi gases avoids some of the complications of lattice 
gauge theory. With nonrelativistic fermions there is no chiral symmetry, so there is no fermion 
doubling problem. Without nonabelian gauge fields, there are no topological sectors to worry about 
sampling with the correct distribution. This is an ideal system for concentrating on Monte Carlo 
methods for fermions. In addition to both being strongly coupled theories, QCD and unitary Fermi 
gas both exhibit spontaneous symmetry breaking at low temperatures, restored at some critical 
temperature. Like finite density QCD, the fermion matrix of the unitary gas can have a nonpositive 
determinant. 

2. Setup 

We start from the Fermi-Hubbard model, which in the grand canonical ensemble reads, 

H = Ho + H l = Y,(£k-l^a)cl a c ka + U^4t c ^ c li c ^ C 2 - 1 ) 

k.cr x 

where £t = ^L] =1 (l —coskj) is the discrete dispersion relation, jx a the chemical potential and 
CkCT (ckcr) the time-dependent fermionic creation (annihilation) operator. We use h = ks = 2m = 1. 
The index a € {t, 1} labels the fermionic species. The coupling constant U < corresponding to 
attractive interaction can be tuned so that the scattering length becomes infinite. The corresponding 
value is U = —7.914. We work on a 3D periodic lattice with L? sites. The continuum limit can be 
taken by extrapolation to vanishing filling factor v = (Lff c xcjc: xfT ) — >• 0. 

According to §\§ the partition function for this model can be written as a series of products of 
two matrix determinants built of free finite-temperature Green's functions. If jJ.^ = jJ.^ (the balanced 
case) these determinants are equal so that all terms in the series are positive and it can be used a 
probability distribution for Monte Carlo sampling. The order parameter for the phase transition is 
given by a two-point correlation function of the operator c x |C x ^. To obtain T c from the numerical 
data, previous work [^, |3|] used a procedure involving an approximation which introduced a sys- 
tematic error. We have improved the data analysis method so that this approximation is no longer 
necessary Our final result will be the dimensionless quantity T c /ep, where the Fermi energy 
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Figure 1: Schematic plot of the average sign near the critical point. The shaded area covers the range of 
values the sign can take at different values of lattice size and chemical potential. The lower boundary of this 
area is the "worst-case" curve of the sign, corresponding to lowest densities and largest lattice sizes used. 

Ep = (37T 2 v) 2 / 3 is the only energy scale of the system. The continuum limit is taken by obtaining 
T c /ep for different values of v and extrapolating to vanishing filling factor [§, @]. 

A detailed description of our numerical setup is given in [Q] and [^J. We use the DDMC 
algorithm as introduced in [Q] with a few modifications which increase the efficiency by reducing 
autocorrelation effects that are present in the original setup. Also we generalise the algorithm to 
the spin-imbalanced case. If jUf ^ jUi, a sign problem arises, since the distribution function is no 
longer positive for all configurations. To deal with this problem we make use of the "sign quenched 
method", which is based on the "phase quenched method" known from lattice QCD [^]. The idea 
is to write the thermal distribution as a product of its modulus and its sign, and to use the positive 
function given by the modulus as the new probability distribution. This implies a reweighting of 
each MC estimator, in this case simply a multiplication with the relative sign of the two matrix 
determinants. Additionally, each expectation value must be divided by the expectation value of the 
sign. If the latter is is close to zero, numerical errors will be very large, as it happens in QCD. 
However, for the unitary Fermi gas the sign remains very close to unity for small imbalances, as 
shown in Fig. [IJ so that sign quenching is applicable for imbalances up to approximately Ajx = 
0.2ep. Our method can provide a useful tool to examine the trend of the critical temperature for 
small deviations from the balanced limit. 



3. Results 

We obtained data at 25 different values (jU-|-,jU .), of which 8 were at ^ = The lattice 
sizes varied between 4 3 for the highest filling factor and 26 3 for the lowest, so that the volume 
range in physical units was approximately constant. As discussed in for the balanced case, the 
dimensionless physical observables scale linearly with v 1 / 3 for sufficiently small v. With our data 
this behaviour is seen for v 1 / 3 ;$ 0.75. This condition was fulfilled for 23 out of the 25 points and 
in particular for 7 out of the 8 balanced points. 

Since the chemical potential difference is less prone to numerical errors, we use A^/sp = \fu — 
jUi |/ Bp instead of the relative density difference Av/v to quantify imbalance. For the values of 
imbalance considered in our study these two quantities are proportional to each other, with Av/v = 
0.l22(2)Afi/£p, see Every physical observable X is a function of filling factor v and imbalance 
h = An/sp. We fit our data to a three dimensional surface, where the following assumptions 
are made for the form of the fitted function: 
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• At fixed imbalance, X is a linear function of v 1 / 3 with slope a^ x \h): X(y,h) = X(h) + 
a^ x \h)v 1 ^ . This is a generalisation of the relation valid in the balanced case. 

• X(h) and a^ x \h) viewed as functions of the imbalance h can be Taylor expanded. 

• Due to symmetry in h all odd powers in the expansions of X(h) and aS x > (h) have to vanish. 
If we expand X(h) and (h) to leading order in h the fitted function becomes 

X(v,h) = X + X 2 h 2 + (af } + of ] h 2 )v l/3 . (3.1) 

In the following we will present our fit results for several physical quantities. The results from 
sections 



3. 1 , 3.2 and 3.3 were presented in more detail in 



3.1 Critical Temperature 

Fitting a line through the points with jitf = jU; results in T c /e F = 0.173(6) - 0.16(l)v 1 / 3 with 
% 2 /d.o.f = 0.39, as shown in Fig. ^[ For comparison we also fit a quadratic through all 8 data 
points, resulting in a continuum value of T c /sp = 0.188(15), which is in excellent agreement with 
the linear extrapolation. This confirms that sub-leading corrections proportional to v 2 / 3 can indeed 
be neglected for sufficiently small v. 

Now we also include data with jx^ ^ jx^. The best fit according to d3~l] ) yields 7b = 0.171(5), 
a^ T) = -0.154(9), T 2 = 0.4 ±0.9 and a ( 2 T) = -0.7 ± 1.9 in units of e F , with # 2 /d.o.f.= 0.43. 
Note that the T 2 value corresponding to the minimal % 2 is positive, which is forbidden by physical 
arguments - the critical temperature can only decrease with increasing imbalance. Since the % 2 
function is very flat along the T 2 direction, forcing T 2 = results in ^ 2 /d.o.f.= 0.44. From the 

(T) 

error on T 2 we derive the lower bound T 2 > —0.5. The best fit values for 7b and a () are in 
excellent agreement with the ones obtained from the fit of the balanced data only. 

(T) 

By simplifying the fit model setting a 2 =0, the lower bound on T 2 can be tightened to 

(T) 

T 2 > —0.04. The other parameters Tq and a$ agree with the results from the previous fit. This 
fit has # 2 /d.o.f.= 0.41 and is still consistent with T 2 = 0. For this reason we also perform a fit to 
constant T c {h) and a^ T \h), and obtain T c (v,h) = 0.1720(45) -0.156(8)v 1 / 3 with ^ 2 /d.o.f.= 0.41. 
Again the result agrees with the previous fits. We also performed fits using the jackknife method 
and several robust fits and obtained consistent results. A three dimensional plot of the data together 
with the constant surface fit is presented in Fig. ^ (left). 

3.2 Energy per particle 

The total energy is composed of the kinetic energy = — ^L XifT Cx CT V 2 c XCT ^) and the inter- 
action energy E- mt = (Hi). For the explicit MC estimators see [Q]. Our results are obtained at T c , 
but the temperature dependence of the energy per particle was found to be weak. Using only bal- 
anced data we obtain the continuum value E/NSp = 0.276(14), or E/Epc = 0.46(2), in units of 
the ground state energy of the free gas Ep q = (3/5)A^£ F . The goodness of fit is £ 2 /d.o.f. =2.1. 
Since with increasing imbalance interactions become suppressed, the absolute value of the interac- 
tion energy must decrease. This in turn means an increase of the total energy, since the interaction 
energy is negative. As we did for the critical temperature we fit the energy in units of Epc to the 

,(E) 



function (gTJ) and obtain the best fit parameters E = 0.440(15), af> = -0.17(3), E 2 = 3.4±2.2 



and = —3.1 ±4.5, with # 2 /d.o.f.= 2.8. These results are consistent with the balanced fit. 
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Figure 2: T c versus v 1 / 3 for different values of 
the chemical potential at ju-j- = jUi . The solid line 
is the linear extrapolation of the 7 data points with 
V 1 / 3 J$ 0.75 (filled circles). The dashed line cor- 
responds to a quadratic fit through all data points. 




T C /£ F 



Figure 3: Projection of the data for the average 
chemical potential onto the (v 1//3 -jti) plane. Red 
circles denote the balanced data and blue triangles 
data at non-zero imbalance. The solid line is the 
constant fit; dashed lines indicate uncertainty. 
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Figure 4: T c /ef with the constant fit (left) and E/Efg with the quadratic fit (right) versus v 1 / 3 and A[i/ef- 

Forcing af ] = yields a best fit result of E(v,h) = 0.444(13) + 1.9(3)/j 2 - 0.18(2)v 1/3 with 
X 2 /d.o.f.= 2.7, which agrees with the previous result. For a plot of the data see Fig. |4] (right). 

3.3 Chemical potential 

For the chemical potential at T c we obtain the continuum value h/sf = 0.429(9) with # 2 /d.o.f. 
= 2.8 using only balanced data. A similar analysis can be performed for the average chemical 
potential \ijep = \H\ + Mj,|/2£f in presence of an imbalance. Since this quantity is not expected 
to depend on the imbalance we fit our data to a constant function and obtain jj,(v,h) = 0.429(7) — 
0.27(l)v 1 / 3 in units of Sp with # 2 /d.o.f.= 1.1. This is in very good agreement with our balanced 
result. A plot of the data and the fit are shown in Fig. |5[ 

3.4 Contact density 

Another important quantity is the contact density, which can be interpreted as a measure of the 
local pair density [j8|]. The contact plays an important role for several universal relations derived by 
Tan [^]. The definition of the contact is C = m 2 goEi nt , where go is the physical coupling constant 
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Figure 5: Left: The contact density versus v 1 / 3 for balanced data together with the linear fit. Right: The 
contact density versus filling factor and imbalance. The surface corresponds to the quadratic fit. 



[@, 10], and it is related to the contact density ^ via C = f ^(r)d ' r, or for homogeneous systems 
simply C = "^V. The dimensionless quantity 'tf /£f can be expressed as 

V/e f = (UE int )/(4L 3 e 2 F ). (3.2) 

Using only balanced data the best fit is -jf/e| = 0.1102(11) - 0.033(2)v 1 / 3 with x 2 /d.o.f. = 1.8. 
In the presence of an imbalance we expect the interaction energy and hence the contact density to 



decrease. A three dimensional fit of the data to (gTTJ) yields C = 0.1101(9), C^ C) = -0.033(2), 
C 2 = -0.15(16) and c^ C) = -0.29(36), with ^ 2 /d.o.f.= 1.5. This is consistent with the bal- 
anced fit. The parameter C 2 is negative as required. Forcing a\ = yields Co = 0.1099(8), 
a { p = -0.0322(15) and C 2 = -0.01(2), with % 2 /&.o.f.= 1.5. Fitting to a constant function yields 
tf(v,h) = 0.1097(8) - 0.0320(14)v 1 / 3 with x 2 /d.o.f.= 1.4. Figure | summarises our results. 

3.5 Comparison with the literature 

Our final result for T c using both balanced and imbalanced data, T c /£p = 0.171(5), is signif- 



icantly higher than the previous result from [g], where T c /Ep = 0.152(7). In [11] the result of [g] 
was found to agree with a continuous space-time DDMC method. The authors of ^ found an 
upper bound of T c /ef is 0.15(1). They used an auxiliary field Monte Carlo approach and extracted 
T c using the same approximation as [||] and [JTTJ] , which might explain the discrepancy between 
our results. Through extrapolating Monte Carlo results of low-density neutron matter, the authors 



of [|12J found T c /Ep =0.189(12), which agrees with our result. There are also results obtained 



with the Restricted Path Integral Monte Carlo method [|T3|], T c /Ep ~ 0.245, and an upper bound of 



T c /ef < 0.14 obtained with a hybrid Monte Carlo method Q14|]. Results from an £ -expansion are 



also available [J15J. For comparison, the critical temperature in the BEC limit is Tbec = 0.218£f . 

Our result for the energy per particle E/Efg = 0.440(15) shows excellent agreement with the 
value E/Efg = 0.45(1) at T c quoted in The value quoted in [@] isE/N£ F = 0.31(1), which 
roughly corresponds to E/Efg = 0.52(2). Our result for the chemical potential \ljEp = 0.429(7) 
differs from \ljEp = 0.493(14) quoted in [Q], but is consistent with the value \xjEp = 0.43(1) 
quoted in Some theoretical predictions for the contact density are also available, but to our 
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knowledge only at temperatures much lower than T c (see |Q] and references therein). 

There are several recent experimental studies of the homogeneous unitary Fermi gas. The di- 
rect measurement presented in [|TJ] is T c /£p = 0.157(15), which agrees well with our result, and 
jJ./Sf = 0.49(2) at T c , which differs from our value. The values from Jl7|], T c /Ep = 0.17(1) and 
jl/SF = 0.43(1) at T c , show excellent agreement with our results. Their result for the energy per 



particle E/Nsf = 0.34(2) at T c is higher than our value. In another experimental work [18] an 
estimate for T c at zero imbalance is extrapolated from data at higher values of imbalance. An ex- 
perimental value for the contact density of the homogeneous unitary Fermi gas at zero temperature, 
tf/ej = 0.1184(64), is presented in 0]. 

4. Outlook and Acknowledgements 

We have presented a Monte Carlo calculation of several thermodynamic observables of the 
unitary Fermi gas with equal and unequal chemical potentials in the two spin components. The 
improved DDMC algorithm with sign quenching also offers the intriguing possibility to explore 
the case of unequal masses of the two species. If ^ m ^ the dispersion relations are different for 
the two components and the mass ratio enters as a new parameter. As in the spin-imbalanced case 
the two matrix determinants no longer need to be identical, so that sign quenching is required. The 
mass ratio is expected to influence the behaviour of the system significantly, so that exploring the 
phase diagram promises many new interesting insights. 

This work used resources provided by the Cambridge High Performance Computing Facility. 
OG is supported by the German Academic Exchange Service (DAAD), the EPSRC and the CET. 
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